## Code that creates Figure 1 in the main text. 
## Last changed: 2025-05-23

sink("C:/Userdata/Shared/Dofiles/DoAnalysis/JanKalle/ReplicationFiles/AnalysisScripts/Output/Figure1.txt", split=TRUE)

library(PSweight)
library(cobalt)
library(MatchIt)
library(distances)
library(moreparty)
library(party)
library(caret)
library(bonsai)
library(randomForest)
library(tidymodels)
library(tidyverse)
library(rmarkdown)
library(tinytex)
library(tibble)
library(knitr)
library(car)
library(fst)
library(stargazer)
library(pryr)
library(stringr)
library(fst)
library(lubridate)
library(data.table)
library(summarytools)
library(ggplot2)
library(cowplot)
library(ggpubr)
library(haven)
library(broom)
library(fixest)
library(modelsummary)
library(gt)
library(MASS)
library(survey)
library(collapse)
library(readxl)
library(Hmisc)
library(DescTools)
library(ranger)
library(stablelearner)
library(future.apply)
library(lme4)
library(partykit)
library(ggiplot)

  gc()
# Merge on earnings data for the years 1982 to 1991 from the Income and Tax register.

# Read in data on the consumer price index to calculate earnings in 1991 years prices.
  KPI <- read_fst("D:/Data/ExtData/KPIData/SCBKPI_5523.fst", as.data.table=T)
  kpi91 <- KPI[Year==1991, MeanKPI]
  KPI[, MeanKPI_91 := MeanKPI/kpi91]
  setnames(KPI, "Year", "Ar")

  IOT <- data.table()
  
  for(i in 1982:2021) {
    iot <- read_fst(paste0("D:/SCB_ConPol/Fst/IoT/IoT_", i, ".fst"), as.data.table=T)
    iot <- iot[, cID := .N, by = LopNr][cID == 1 & !is.na(LopNr), ][, cID := NULL]
    iot[, DAkassa := as.numeric(TKASSA>0)]
    iot[, Ar := i]
    iot[, ArbInk := CARB]
    iot <- iot[, .(LopNr, Ar, DAkassa, TKASSA, ArbInk)]
    
    if(i>1993) iot[, ArbInk := ArbInk/100]
    iot[, ArbInk := ArbInk/10]
    print(qsu(iot))
    IOT <- rbindlist(list(IOT, iot), use.names=T, fill = T)
    print(i)
  }
  
# Specify the name of the treatment variable
  tvar <- "T94_90"

# Specify the the birth cohorts to be used for the analysis
  lbyr <- 1942
  ubyr <- 1967

# Specify if only include Swedish citizens 
  citizen <- 1

# Specify the propensity score specification
  psformula <- formula(T94 ~
                         Kon + DAStod_1982 + DAStod_1983 + DAStod_1984 +
                         DAStod_1985 + DAStod_1986 + DAStod_1987 + DAStod_1988 +
                         DAStod_1989 +
                         DAntBarn_91 + Sambo_91 + InvBak  |
                         SSYK3_90 + SNI2_91 + Kommun_91 +
                         rCARB_1982 + rCARB_1983 + rCARB_1984 + rCARB_1985 +
                         rCARB_1986 + rCARB_1987 + rCARB_1988 + rCARB_1989 +
                         pArbInk_1990 + pArbInk_1991 + UtbAr_91 +
                         FodAr + Nom82 + Nom85 + Nom88 + Nom91)
  
  
  NomVald <- read_fst("E:/ProjData/Jan_Kalle/JoPReplications/DB_NomVald.fst", as.data.table=T)
  db <- read_fst("E:/ProjData/Jan_Kalle/JoPReplications/DB_Nom.fst", as.data.table=T)
  
  db[, T94 := get(tvar)]
  mdb <- db[between(FodAr, lbyr, ubyr) & (SvMedb82_18 >= citizen) & AllaVal == 1,]
  rm(db)
  gc()
  
  # Use logit model to estimate propensity scores
  mdb <- mdb[complete.cases(mdb),]
  tm <- proc.time()
  logmod <- feglm(psformula,  
                  data = mdb, family="binomial")
  proc.time()-tm
  summary(logmod)
  
  # Extract predicted probability (propensity scores)
  if(logmod$nobs == logmod$nobs_origin){ 
    mdb[, pT94 := predict(logmod, type = "response")]
  }else{
    mdb[logmod$obs_selection$obsRemoved, pT94 := predict(logmod, type = "response")]
  }
  
  mdb <- mdb[complete.cases(mdb),]
  
  # Now use the Full matching approach of Sävje et al. to obtain ATT and ATE weights.  
  
  tmp <- proc.time()
  set.seed(7892)
  m.out <- matchit(T94~1, data = mdb,
                   distance = mdb$pT94, method = "quick", estimand = "ATT")
  proc.time()-tmp
  m.data1 <- as.data.table(match.data(m.out))[, .(LopNr, weights)]
  setnames(m.data1, "weights", "gm_att")
  
  # Extract the matching weights and merge them onto the main data set. 
  mdb <- m.data1[, .(LopNr, gm_att)][mdb, on = "LopNr"]
  gc()
  IOT <- mdb[, .(LopNr, gm_att, FodAr, T94)][IOT, on = "LopNr"]
  IOT <- IOT[!is.na(gm_att),]
  rm(mdb)
  gc()
  
  # Calculate earnings in fixed prices (1991)
  IOT <- KPI[, .(Ar, MeanKPI_91)][IOT, on = "Ar"]
  IOT[, fp_ArbInk := ArbInk/MeanKPI_91]
  IOT[, Alder := Ar-FodAr]

# Exclude individuals over 65 and calculate variable means for the various outcome variables using matching weights.  
  mIOT <- IOT[Alder<66,]
  
  # Labor earnings in fixed 1991 prices by year and treatment status
  mIOT[, wLon := qsu(fp_ArbInk, w=gm_att)[2], by = .(T94, Ar)]
  
  # Share of individuals receiving unemployment benefits by year and treatment status.
  mIOT[, wArbLos := qsu(DAkassa, w=gm_att)[2], by = .(T94, Ar)]
  mIOT[, occ := seq_len(.N), by = .(T94, Ar)]
  
# Plot the data and save the results for Figure 1.   
  g0 <- ggplot(mIOT[occ==1,], aes(Ar, wArbLos, linetype=as.factor(T94))) +
    geom_point() + 
    geom_line() +
    scale_x_continuous(breaks = seq(1980, 2022, 5),  labels = seq(1980, 2022, 5), limits = c(1980, 2022)) +
    scale_y_continuous(limits = c(0, 1 ), breaks = seq(0, 1, by = .1)) +
    xlab("Year") +
    ylab("Unemployment Benefit Incidence") +
    theme_classic(base_size = 18) +
    scale_linetype_manual(labels = c("Non-Job Losers", "Job Losers"), values = c(1, 2), name = "") +
    theme(legend.position = "bottom", panel.grid.major.y = element_line(color = "grey95"))
  g0
  
 g1 <- ggplot(mIOT[occ==1,], aes(Ar, wLon, linetype=as.factor(T94))) +
   geom_point() + 
   geom_line() +
   scale_x_continuous(breaks = seq(1980, 2025, 5),  labels = seq(1980, 2025, 5), limits = c(1980, 2025)) +
   scale_y_continuous(limits = c(0, 300), breaks = seq(0, 300, by = 50)) +
   xlab("Year") +
   ylab("Yearly Earnings (1000 SEK)") +
   theme_classic(base_size = 18) +
   scale_linetype_manual(labels = c("Non-Job Losers", "Job Losers"), values = c(1, 2), name = "") +
   theme(legend.position = "bottom", panel.grid.major.y = element_line(color = "grey95"))
 g1
 
 ggarrange(g0, g1)
 ggsave("C:/Userdata/Shared/Dofiles/DoAnalysis/JanKalle/ReplicationFiles/AnalysisScripts/Figures/Figure1.pdf", width = 15, height =7.5, units = "in") 
 ggsave("C:/Userdata/Shared/Dofiles/DoAnalysis/JanKalle/ReplicationFiles/AnalysisScripts/Figures/Figure1.eps", width = 15, height =7.5, units = "in") 
 
sink()